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Different computational scheme for calculating surface integrals in anisotropic Brillouin zones 
are compared. The example of the transport distribution function (plasma frequency) of the ther- 
moelectric Material Bi2Te3 near the band edges will be discussed. The layered structure of the 
material together with the rhombohedral symmetry causes a strong anisotropy of the transport 
distribution function for the directions in the basal (in-plane) and perpendicular to the basal plane 
(out-of-plane) . It is shown that a thorough reciprocal space integration is necessary to reproduce 
the in-plane/out-of-plane anisotropy. A quantitative comparison can be made at the band edges, 
where the transport anisotropy is given in terms of the anisotropic mass tensor. 

PACS numbers: 71.18.+y,71.20.Gj,72.20.Pa,72.80. Jc 



I. INTRODUCTION 

Thermoelectric materials are nowadays 
widely used to convert waste heat into electri- 
cal energy or for cooling purposes. 1 The main 
advantage of this conversion process is the 
absence of moving parts or liquid or gaseous 
components. The efficiency of the process is 
rather low and is mainly determined by the 
dimensionless parameter 



ZT = 



aS 2 



called figure of merit. It is determined by the 
specific conductivity a, the thermopower 5, the 
absolute temperature T, and the thermal con- 
ductivity ft, comprising the electronic and lat- 
tice parts K e and kl- 

New experimental techniques allow for 
the preparation of nanostructered and low- 
dimensional thermoelectric devices which are 
supposed to possess larger ZT values Val- 
ues up to 2.4 are reported for multilayered 
Bi2Te3/Sb2Te3 systems^ A microscopic un- 
derstanding of these effects can be obtained by 
calculating the conductivity and powerfactor in 
the semi-classical limit exploiting the relaxation 
time approximation. The so-called transport 
distribution function has to be determined as 
a function of electron energy E as integrals 
of surfaces of constant electron energy in re- 
ciprocal space.— If one considers an anisotropic 
material, like Bi2Te3 , these quantities show a 
strong directional dependence, which can sup- 
port the enhancement of figure of merit in ther- 
moelectric heterostructures. The calculation of 
these isoenergetic surface integrals requires a 



thorough integration in k space. The impli- 
cations of different integration scheme will be 
demonstrated. To this end the paper is orga- 
nized as follows. We start with a derivation of 
the transport anisotropy near the band edges 
as a function of the effective mass tensor. This 
provides benchmark numbers at certain ener- 
gies to check the numerical integration schemes. 
The systematic deviations of the integration 
schemes are related to the anisotropic structure 
of the reciprocal space of the rhombohedral lat- 
tice and the anisotropy of the band edges con- 
cerning the effective masses are discussed af- 
terwards. At the end, a simplified analytical 
model for a two-dimensional integration will be 
discussed to show the influence of the rhombo- 
hedral lattice anisotropy on the transport quan- 
tities. Results for a free electron model show 
a strong dependence on the procedure to fill 
the reciprocal space with tetrahedrons in a way 
that the symmetry operations of the lattice can 
still be applied to reduce the numerical effort. 



II. INVERSE MASS TENSOR AND 
CONDUCTIVITY ANISOTROPY AT 
BAND EDGES 



We determined the band dispersion for 
Bi2Te3 at the experimental lattice con- 
stants with the atomic positions taken from 
literature^ The topology of the band structure 
is described elsewhere. 8 Assuming a constant 
relaxation time the energy dependent matrix 
valued transport distribution is defined as£ 



(27r) 3 ft 



e(k)=E 



dS 



v a (k)vp(k) 



(1) 

with a, j3 the Cartesian coordinates, k a com- 
bined index of reciprocal space vector k and 
band index e{k) the band energy, and v a (k) 
the group velocity in the direction a. To ob- 
tain the anisotropy cf xx /cf zz as a point of ref- 
erence, we assume in the vicinity of the band 
edges parabolic bands in an anisotropic effec- 
tive mass model. The inverse mass tensor M 
close to the band edges is diagonalized in the 
form 



band structure on a very dense mesh in k-space 
corresponding to 400 points along a reciprocal 
lattice vector to obtain convergence concern- 
ing the position of the extremum and the in- 
verse mass tensor. The values for the valence 
and the conduction band are summarized in 
Tabs.Il and IP 



VBM 


M % 






Ci,z 


1 


-41.3 


0.5000 


-0.8660 


0.0000 


2 


-7.42 


0.6000 


0.3463 


0.7212 


3 


-0.507 


0.6246 


0.3606 


-0.6927 



M 



m 



diag(Mi,M 2 , 



^i^y^y 



M 3 ), 
z (i 



TABLE I: Inverse effective mass tensor eigenvalues 
Mi and eigenvectors c% at the valence band maxi- 
mum (VBM). 
with eigenvectors 



1...3) 



with e a the basis vectors of the Cartesian co- 
ordinate system. The transport distribution for 
the anisotropic effective mass model along the 
main axes q of the effective mass ellipsoid are 
proportional to 



an oc 



m 



m 



oc — = rriM; with 



(3) 



(4) 



where the masses have to be chosen positive 
for both the valence band maximum and the 
conduction band minimum. 

The following expressions for the in-plane 
and cross-plane conductivities project the mass 
tensor of general orientation given by eq. ([2]) to 
the Cartesian coordinates. 



oc 



oc 



oc 



Me, 



(ci )X ) 2 Mi + (c 2 , x ) 2 M 2 + (c 3 , x ) 2 M 3 
) 2 M 1 + (c 2 ,y) 2 M 2 + (c 3 , v ) 2 M 3 
(c hz ) 2 M 1 + (c 2>2 ) 2 M 2 + (c 3)2 ) 2 M 3 



(5) 



Due to the space group D% d (R3m) the con- 
sidered band extrema are two- and six-fold de- 
generate. In the case of the rhombohedral lat- 
tice this summation leads to equal contribu- 
tions a xx and a yy and keeps the a zz unchanged. 
The in-plane transport distribution is given for 
symmetry reasons by cry = (a xx + o yy ) /2. and 
the cross-plane component by a± = a zz . In 
the following the term <j xx will be used synony- 
mously for the in-plane component cry, so the 
anisotropy is denoted as cf xx /cf zz . 

The mass tensor in Bi2Te3 is parameterized 
close to the band edges using the calculated 



(2) 



The position qs of the six-fold degenerate va- 
lence band maximum in units of inverse Bohr 
radii is (0.372199,0.644655,-0.0299675) on the 
plane (rZU). The effective mass ellipsoid is 
very anisotropic. The angle (j) of the long axis 
of the ellipsoid with the (xy) basal plane is 
43.8°, which is in good agreement to other cal- 
culations and experiments, which show a quite 
spread^ - — The transport anisotropy ratio de- 
termined by the effective mass tensor for hole 
states close to the valence band maximum ac- 
counts to (T xx ja zz = 5.40. 
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5.62 


1. 
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5.62 
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1. 


0. 
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1.20 
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0. 


1. 



TABLE II: Inverse effective mass tensor eigenval- 
ues Mi and eigenvectors Ci at the Conduction band 
minimum (CBM). 

Tabs. Ull summarizes the effective mass ten- 
sor at the two-fold degenerate conduction band 
minimum with a the position qs in units 
of inverse Bohr radii (0,0,0.0552), which is 
about one third of the line(rZ). The transport 
anisotropy is cr xx /a zz = 4.67. 

With these values the convergence of the 
transport distribution given by eq. (PQ) deter- 
mined by interpolation schemes in the Brillouin 
zone for energies close to the band edges can be 
quantified. 

Using the calculated band structures and the 
derived effective mass tensors we obtain an 
anisotropy cr xx /a zz of 5.40 and 4.67 at the va- 
lence and conduction band edge, respectively. 
These values are marked by dots in Figs. [T]) 
and[2j). 



III. CONDUCTIVITY ANISOTROPY: 
COMPARISON OF DIFFERENT 
INTERPOLATION SCHEME 

The transport distribution cr(E) of Bi2Te3 is 
calculated by two methods. The main distinc- 
tion is the determination of the group velocities 
v(k). 

The tetrahedron method divides the irre- 
ducible part of the Brillouin Zone (BZ) into dis- 
joint tetrahedra. The group velocity is obtained 
by a linear interpolation of the band energies 
at the 4 corner points and approximates v(k) 
in the volume of the tetrahedron. 12 The sec- 
ond method determines the velocities as deriva- 
tives along the lines of the Blochl mesh^ in the 
whole Brillouin zone. The directions of these 
lines are parallel to the reciprocal space vectors 
and so the anisotropy of the real lattice is re- 
flected in these vectors. For Bi2Te3 with a large 
ratio c hex /a hex of 6.95 7 the real space unit cell 
is very prolongated and the reciprocal lattice 
vectors are quite close to the (xy) basal plane 
with an angle 6' — 118.05° between them very 
close to the maximum value of 120°. Project- 
ing back these so-called mesh velocities to the 
Cartesian components quite large errors occur 
in the resulting velocities as discussed below. 
The results of the integration of eq. (pQ) using 
both interpolations schemes and different den- 
sities of k-points are compared in Figs. [I]) and 
0). 

We start the discussion with the valence band 
maximum because the principal axes of the 
mass tensor are directed along the x, y and z 
axis. The left hand panel of Fig. [1] shows the 
parabolic dispersion along the three principal 
axes with the very different effective masses, 
spreading by a factor of about 80. The right 
hand panel summarizes the anisotropy ratios 
o'xx/o'zz as a function of energy. It is obvious 
that both interpolation schemes give system- 
atic deviations from the expected values. For 
very large k-mesh densities defined by 384 mesh 
points along a reciprocal lattice vector the re- 
sults converge to the correct value at the band 
edge. The necessary densities of k-points are 
very demanding for realistic band structure cal- 
culations with some atoms in the unit cell. The 
tetrahedron method overestimates systemati- 
cally the anisotropy ratio, whereas the mesh ve- 
locity method underestimates the values. This 
behavior will be discussed below for a simpler 
two-dimensional lattice and the found trends 
are confirmed. 

A similar behavior is obtained for the conduc- 
tion band minimum as shown in Fig.[2j). Here, a 
second challenge appears with a local maximum 



of the conduction band at the T point very close 
in energy to the conduction band minimum. In 
addition, a large transport anisotropy is ob- 
tained due to the occurrence of saddle points 
in the band structure, as discussed elsewhere. 8 
The obtained anisotropics for the trans- 
port distribution can be strongly influenced 
by the k-mesh density especially close to the 
band edges. Taking the converged transport 
distributions T,(E) we obtained for small p- 
doping an anisotropy of about 6, and for 
small n-doping an anisotropy of about 9. 14 
This about a factor of 2 larger than found in 
experiment 15 and other calculations^ - — As- 
suming an anisotropy of 0.47 for the averaged 
relaxation time of states traveling along and 
across the basal plane we obtain a good agree- 
ment with experiment. The reason for this scat- 
tering anisotropy has to be elucidated by calcu- 
lations of the microscopic transition probabili- 
ties caused by defects or other scattering cen- 
ters. 



IV. INTERPOLATION SCHEMES: 
TETRAHEDRON VS. MESH 
VELOCITIES 

A. Mesh velocity method 

In the following the implications of the differ- 
ent interpolation schemes for the group velocity 
are discussed. We start with the discussion of 
the so-called mesh velocities. They are deter- 
mined by a numerical derivative of the band 
dispersion along the directions of the recipro- 
cal lattice vectors which span the Blochl mesh. 
The Cartesian components of the velocities are 
obtained by a non-orthogonal transformation. 
If the angles between the basis vectors are very 
large (the maximum is 120°) small errors can 
be largely enhanced. This appears especially 
at band structures with a strong anisotropy 
between directions in (xy)-plane and the z- 
direction. 

Fig. [3] shows schematically the band disper- 
sion of Bi2Te3 close to the conduction band 
minimum by isoenergetic lines. The Brillouin 
zone is shown in reduced size but the angles 
correspond to the considered case of Bi2Te3 . 
As it is obvious, the interpolation of the veloc- 
ities deviates strongly from the correct values 
for k-points close to the band extremum. The 
directions of the reciprocal basis vector mainly 
scan along the (xy)-plane. If the anisotropy of 
the band dispersion e(k) is very strong as in 
the case considered, the mesh velocities tend 
to equalize the in-plane and out-of-plane com- 




FIG. 1: Band structure of Bi2Te3 (left) and conductivity ratio (right) near valence band maximum, group 
velocity interpolation by tetrahedron method (tet) and mesh velocity concept (mesh), number of k-points 
along a reciprocal lattice vector is given as parameter. Observe that the vertical axis gives the energy and 
the anisotropy a xx /a zz is on the horizontal axis. 



ponents of the velocity v XiV « v z which leads 
to an anisotropy closer to unity. This can be 
seen in Figs. [1]) and [2j) by the curves labeled 
'mesh'. This effect is most pronounced close to 
the band edges. 

A quantitative comparison of the accuracy 
of the mesh velocities in given in fig. HI To 
compare the mesh velocities with the exact val- 
ues we have chosen the TZ line and consider 
the z-component of the velocity v z (k z ). The 
minimum of the energy dispersion is the global 
minimum of the conduction band. A local max- 
imum appears at the T point. For sparse k- 
mesh densities a strong deviation of the mesh 
velocities have to be stated. To reproduce the 
point of inflection of the band dispersion by the 
sign change in the velocity a very dense mesh is 
necessary. The velocity in z-direction is overes- 
timated by this method, if the velocities in the 
(xy) directions are larger than in z-direction, 
as illustrated in fig. [3] and typical for Bi2Te3 
in both the valence and conduction band. As a 
result the transport anisotropy cr xx /a zz , which 
is de facto the ratio of the velocities squared, is 
shifted to unity. An increase of the anisotropy 
towards 1 is expected for the opposite case with 



larger velocities along z-direction than along 
(xy) directions. 



B. Tetrahedron method 

The capability of the tetrahedron method to 
calculate the anisotropy of the transport distri- 
bution will be evaluated for a two-dimensional 
mesh for the cases of linear and quadratic inter- 
polation of the velocities. The results for a free 
electron dispersion with two different arrange- 
ments of tetrahedrons in the irreducible part of 
a three-dimensional meshes will be compared 
with. 

The regular filling of the two-dimensional lat- 
tice is characterized by the angle a between the 
border line of the triangles Si and S2 and the 
basal (xy)-plane, see Fig. [5] Depending on 
the interpolation method the velocities have to 
be determined for the triangles Si and S2 with 
a linear scheme and can be determined at the 
points A, B, and C with a quadratic interpo- 
lation scheme. The second task requires the 
knowledge of the energy dispersion e(k) at ad- 
ditional points on the edges. Solving the lin- 
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FIG. 2: Band structure of Bi2Te-3 (left) and conductivity ratio (right) near conduction band minimum, see 

fig. El 



ek = const. 




FIG. 3: Determination of the group velocity by 
means of so-called mesh velocities vm along the 
reciprocal basis vectors. Anisotropic band dis- 
persion is sketched by elliptical isoenergetic lines 
€k = const, in the (k x ,k z )-p\axie. The rhombohe- 
dral Brillouin zone for the case of Bi2Te3 is shown 
by thin lines. 



ear set of equations for the Cartesian compo- 
nents of the velocities v\ and V2 for a free elec- 
tron dispersion, one finds both velocities equal 
pointing in the same direction along the bor- 
derline of the two triangles v\ \ \v2\ | (cos a, sin a). 
So, the anisotropy of the transport distribution 
&xx/&zz is just given by the ratio of the com- 
ponents v x and v z squared 



Cxx/ O ' zz — 



cos 2 a 



sin 2 a tan 2 a 



(6) 



For the given c/a ratio of Bi 2 Te3 the angle a 
is about 23° which results an anisotropy value 
of 5.8 which is much to large in comparison of 
the expected unity ratio. Improving the inter- 
polation scheme of the velocity to second order 
the velocities va, vb, and vc are determined 
correctly and the error of the transport distri- 
bution function is mainly determined by the ap- 
proximation of the Fermi surface (a line in the 
two-dimensional case) by the two line segments 
a\ and d2- Due to the canceling of prefactors in 
the ratio <j xx /<j zz it is sufficient to consider the 
sum of the segment lengths \a\\ and |a,2 1 times 
the square of the corresponding velocity com- 
ponents of va and vb for area Si, and vb and 
vc for area S2 , see eq. (pQ) 



a xx ex cos 2 ay (1 — sin 2 a) 2 + cos 2 a 
+ (1 + cos 2 a) J (1 - cos 2 a) 2 -{ 



sin 2 a 



ex cos 2 aV 2 — 2 sin a + (1 + cos 2 a) a/2 — 2 cos a 
ex sin 2 ol\J2 — 2 sin a + (1 + sin 2 a)\/2 — 2 cosa(7) 



Bi 2 Te 3 a aTe =4.384A (rhombo) 
0.344 -i 
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FIG. 4: Band structure e(k z ) of Bi2Te3 on the line 
T-Z (top) and mesh velocities in z direction (bot- 
tom) for different densities of the k-mesh. 'drv' de- 
notes the exact result from the derivate de(k z )/dk z . 
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FIG. 6: Anisotropy of conductivity for a free elec- 
tron model in dependence of the assumed underly- 
ing rhombohedral lattice characterized by the an- 
gle a between the reciprocal lattice vectors and 
the (xy) basal plane. As comparison the conduc- 
tivity anisotropics of the 2D model (linear and 
quadratic interpolation) are shown with solid lines. 
Top panel: k-mesh created along (111) directions, 
bottom panel: k-mesh created along (001) direc- 
tions. 
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FIG. 5: Schematic sketch of an anisotropic k-mesh 
division characterized by the angle a. The con- 
stant energy surface is approximated by the poly- 
gon ABC with the parts a\ and a2 in the triangles 
Si and S2. va, vb, and vc denote the velocities. 



From these expressions the ration cf xx /cf zz 
is evaluated and shown as (thicker) blue lines 
in fig. [5] together with the values obtained by 
the linear velocity interpolation as (thinner) red 
lines. The linear scheme shows strong devi- 



ations especially for very anisotropic lattices 
with small or large angles a. This error is 
much reduced by the second order interpolation 
scheme with a maximum error of about 30%. 

Performing these procedures in a three- 
dimensional lattice requires the set-up of k- 
point mesh filling the irreducible part of the 
Brillouin zone and allowing for a disjunct tetra- 
hedron arrangement. Two filling schemes are 
evaluated, which are based on cubes (in a given 
basis) which are filled with 6 tetrahedron each 
if completely inside the irreducible part, oth- 
erwise they are partially used to fill the irre- 
ducible part. The first one is based on cubes 
with the main axes directed along the (111) 
reciprocal lattice directions. The results for 
the anisotropy cf xx /cf zz with a linear and a 
quadratic velocity interpolation are shown as 
crosses and plus signs in the upper panel of 
fig. [6j respectively. To define a similar angle Z 
characterizing the anisotropy of the rhombohe- 
dral reciprocal lattice as in the two-dimensional 



case we have chosen (ir — Z(gi,#2 + 93))/% to 
ensure the isotropic simple cubic lattice to be 
characterized by an angle 7r/4. 

The second filling scheme uses the (001) di- 
rections as basis. The results are shown in the 
lower panel of fig. [6] The results of the two- 
dimensional model can be partially reproduced, 
especially the case of a isotropic lattice at an an- 
gle of 45°. On average the deviations are larger 
for the linear interpolation scheme in compari- 
son to the quadratic one as expected from the 
two-dimensional model. The partially opposite 
behaviour of the three-dimensional integration 
scheme with respect to the two-dimensional one 
needs further clarifications. 

The results in fig. [6] are evaluated for ener- 
gies very close to the free electron band min- 
imum. In these cases the isoenergetic surface 
is approximated by a few trigonal elements 
only and the largest anisotropics caused by 
the interpolation errors of velocities and sur- 
face areas are expected. For larger energies 
these discrepancies disappear quickly. The dis- 
cussed deficiencies of the interpolation schemes 



in k-space are restricted to energies close to 
band extrema, which appear in transport prop- 
erties of medium-doped semiconductors. In 
cases of very small doping the application of 
an anisotropic effective mass model seems to 
be more advantageous. 



V. CONCLUSIONS 

By means of model and realistic band struc- 
ture calculations for crystals with rhombohe- 
dral symmetry we have shown that the deter- 
mination of the transport distribution function 
requires very dense meshes in k-space. Two 
different methods to determine the group ve- 
locities are evaluated. It is found that they 
underestimate and overestimate the transport 
anisotropy in a systematic manner. For very 
prolongated unit cells the anisotropy in k-space 
requires a thorough check of the convergence 
with respect to k-space density. 
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